For Reference 


= 
© 
=) 
a= 
2 
em) 
<a 
= 
=) 
4 
ey 
Zz 
co) 
3 
= 
<a) 
fel 
ae) 
io) 
= 
= 
fo) 
Zz 

















Ex wars © 
UNIMRSUTATIS 
AUBERUAENSIS 














a ————— 


a : - wT 7) 
7 7 ; 





THE UNIVERSITY OF ALBERTA 


THE APPROXIMATE DISTRIBUTION OF THE REGRESSION COEFFICIENT 


BETWEEN TWO STATIONARY, LINEAR, MARKOV SERIES 


by 


RICHARD ROUTLEDGE 


A THESIS 
SUBMITTED TO THE FACULTY OF GRADUATE STUDIES AND RESEARCH 
IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE 


DEGREE OF MASTER OF SCIENCE 


DEPARTMENT OF MATHEMATICS 


EDMONTON, ALBERTA 


FALL, 1972 





Pus IAI Ts. 


THE UNIVERSITY OF ALBERTA 


FACULTY OF GRADUATE STUDIES AND RESEARCH 


The undersigned certify that they have read, and 
recommend to the Faculty of Graduate Studies and Research, 
for acceptance, a thesis entitled "The Approximate Distribu- 
tion of the Regression Coefficient Between Two Stationary, 
Linear, Markov Series" submitted by Richard Routledge in 
partial fulfilment of the requirements for the degree of 


Master of Science. 


a ——_——— 
eal 


Digitized by the Internet Archive 
In 2021 with funding trom 
University of Alberta Libraries 


https://archive.org/details/Routledge1972 


ABSTRACT 


An approximate density function for the classical 
regression coefficient when the residuals and the "con- 
trolled variable" are independent, Gaussian Markov series 
is derived and is compared to the classical density func- 


tion for the regression coefficient. 
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CHAPTER I 


INTRODUCTION 


HS Introduction and Summary 
The simple autoregressive model, 


ve = Oy So ByX; + Cys 


eagagy ARG 


where x5 


and Si Ahan + Nas 


and the e's and n‘'s are independent, standard normal random 
variables, frequently occurs in the analysis of time series. 
Several methods have been proposed for estimating the para- 
meter, By (see, e.g., Malinvaud (1966)). Suppose, however, 
that the classical, least squares procedure is used; that 
is, suppose that, in the case where E(y;) and E(x,) are 
known, the origin is shifted to make them 0, and By is esti- 


mated by b, = (yy ae ye and in the case where the 


1 
means are unknown, By is estimated by 


b, = E(x, -%) (yg-¥) 1/2 (x,-%)*. 


t All summations, unless explicitly stated otherwise, are 


from i equals 1 to n. 
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What properties do these classical estimators still have? 


Watson (1951) has treated a similar problem, multiple 
linear regression with deterministic explanatory variables. 
He found that the classical estimators were unbiased and 
consistent. Unbiasedness is not hard to show in the model 
under discussion here. 

First note that 


Ex, (6) x; +e;) iz ye estate 


and that 


*  REXiy, — (2x,) (Ey,) 


b. = 
1 2 2 
n=x; - (Ex; ) 


NIX, (a,+8,%; +e; ) = (3x, ) Z(a,+8,x,+e;) 
it eee A. SO ke 
; nix, a (2x, ) 


‘ nix,e; - (Zx,) (Ze, ) 
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iu 
nyx. - (Ix; ) 


Thus cies only contribution By makes to the distributions of 
* , * : 
by and by is a shift, and hence, by and by are unbiased for 
& 
By athe (b, -B,) and (b,-8)), both of which are independent of 
Ba, have zero expectation. Furthermore, since the last two 


statistics are used in testing hypotheses and setting up 


confidence intervals, we can restrict our attention to the 
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case where By = 0. Moreover, both by and by are independent 


of Qa and henceforth we may assume, without loss of genera- 


LeCy co ac a, = 0. With these restrictions, observe that 
Vea S4t 6 


E(b,) = E{E{b,|x,,---.*,}} 


-_ 2 
= E{E{ (2x,y,)/(2x>) |X, 7--- x, F} 


dl 


= E{ (x4) 7 XIX, E(y,)} 


co th 


i”) 
5 
As 
ie] 
oF 
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* 
= E{E{b, |x,,---/*,}} 
= B{E{ntxjy, - (2x,) (2yj)]/loxy - (0%4)4) 
[xy 7-6 7X3} 


Al 


= B{ (nix? - (2x, ) 717 x [nix E(y;) in (Sx, ) (ZE(y;)]} 


= Ops 
Hence, by and De are unbiased. Furthermore, using the fact 
that the x's are independent of the y's, and that the x's 
are symmetrically distributed about zero, intuitively, one 


* 
would expect by and by to be symmetrically distributed about 


zero. 
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In the next section of this chapter, the distri- 
butions of by and b, are derived in the classical case of 
independent observations. Chapter II deals with the more 
general case of non-zero serial correlations. McGregor's 
(1962) derivation of the approximate joint density func- 
tion for the two regression coefficients for the "known 
means" case, and Bielenstein's (1963) corresponding result 
for the "fitted means" case are derived in the first four 
sections of Chapter II. In the fifth section, these joint 
densities are integrated to obtain the approximate densities 
of b, and by. It is then shown that these densities re- 
duce to the classical distributions of section 1.2 when 
eS fy SDs es the sixth section, approximate variances 


are calculated. 


The results of a numerical renormalization with 
some sample graphs are presented in Chapter III, along 
with the results of a simulation to test the accuracy of 


the approximation. 


Programs used in the numerical renormalization 


and simulation are presented in the appendix. 
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1.2 The Classical Densities of the Regression Coefficients 


Suppose a simple random sample of n observations is chosen 


from a bivariate normal population with mean, UL = (Uy Ho), 
2 
oF P05, 95 
and covariance matrix, ‘ 
2 
Bota? oi) 


If the mean vector is known, then it is subtracted from all 
the observations. Denote the resulting sample by 
(x) +z) prees (x +Y,) © 


Then, the joint distribution of the sample is 


ar, a = [210,0, (1-97) 77)" 
EP GLO Y ME CGS 
51 Be RL, Lys 
ae a dx,---dy, 
2(izpe) |: OF 5155 O5 


ete 4 2 1/2 -n 
15 [270,9, (1 Die) ] 


2 2 
. —n Sy 2prs,S5 S5 : 
exp| ———5- | - tt dx,---dy, , (ieee) 
2(l1-o') OF 5455 O5 


2a 2 2am 2) a 
where sj = (Zx5)/n, S505 (Zy<)/n, and r = (Zx-y;)/ns, Sy : 


To obtain the joint distribution of the statistics, Syr Sor 
and r, we need to transform the differential volume element, 
dx,..-.dy, , to an expression depending on Si Sor 4 ds), aso, 
and dr. A geometric approach, first used by Fisher (1915) to 
obtain a corresponding result for the case where u is unknown, 


is most appropriate. 
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Por efixed X,-10 ands , y is con- 
strained to lie on the circle, C, 


of radius V3s sin @. 


(b) 2 








If yr is allowed to vary to r + dr 


.s 


>_ 
ee oe 
ie 
a eRe aa 
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y can lie anywhere on the annular 


(c) section, A, of the surface of the 





sphere with centre O and radius 


¥3s . 
2 


If, furthermore, s is allowed to 
2 


vary to s + ds , Y can lie any- 
2 2 
where in the toroidal region, T. 





VA 


ne 





Figure l: The Differential Volume Element for n= 3. 
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Bincegtre= (Ex, y,)/L (Exe) (fyg) 17/7 | x can be considered 


as the cosine of the angle, say 8, between x and y. Further- 


ona / 2 


more, yns, = [rx] 


= the length of x, and Yns, = the 
length of y. (See Figure 1(a), for the three-dimensional 
case.) 

Piexecnuerearce li Xca, then y can vary over the cone of 
Ang tee omLOBXem(oecer igure UGsprage tbe So is also fixed, then 
y is constrained to range over an (n-2)-dimensional spherical 


2)1/2 


surface of radius Yns.sin(6) = Yns, (1-r If the correla- 


2 
tion is allowed to vary from r to r+tdr, then y can lie any- 
where on a spherical shell of thickness ynsd = Yas, (1-7) */7ar, 
(See Figure l(c).) Finally, if the length of y is allowed to 
vary from yns. cs Yyn(s,+ds.), then the set of allowable 
values of y has a further dimension, say a width, of size ~ ds) 


Hence, for fixed x, 
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dy >. dy, « [seer oes) “SeKier) a drds. 
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(n=s)Y 2 


aed 2 
=ao5 (l-r") ards, 


For fixed Sj, x is constrained to an (n-1)-dimensional 
spherical surface of radius ynSy, and hence of (n-1)-dimen- 


s®-1, tf the length of x is allowed 


sional volume « (Vas, )"*« 

to vary from ynsy eis, Yyn(s,+ds,), x can be anywhere on a 

spherical shell of approximately the same radius and of 
‘ . n=l 

thickness ynds,- Hence, dx, ++ -dx, = 8) ds, 


Nelen=27 2. in=3)/2 
and dx,-..dx dy,---dy, « S$, So (l-r™) ds,dsodr. (U2 20 


Now, using (1.2.1) and (1.2.2), the joint distribution 
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of Syr Sor and r is given by 
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Transforming the variables Si; Sor and r to Sy Sor and 


b, = rso/S), with Jacobian, 
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al 0 0 
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the joint distribution of Sy Sor and by is found to be 
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Noting that r” < 1, and hence, bist /ss i, 


|A 


and SS Sy the yoint distribution of by and Sj) is found 
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ubstituting u = S5-S)by, with du = 2sods,, the integral in 
(1.2.4) becomes 
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Finally, integrating (1.2.7) over the range of Sy: the 


distribution of by is obtained as 
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To evaluate the constant, first make the substitution, 
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noting that the integrand is proportional to the t-density, 
with n degrees of freedom. Hence, (see, e. g., Kendall and 
Stuart (1958), equation (11.44)) 
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CHAPTER II 
APPROXIMATE DENSITIES FOR THE REGRESSION COEFFICIENTS 


2.1 The Form of Two Relevant Generating Functions 

To find the joint densities for the regression coeffi- 
cients, and hence the marginal densities, the joint density 
of the x's and y's is first calculated. From this result, 
approximate expressions for the joint moment generating func- 
tions of Ba Ex;Yqe and by? for the ‘known means' case, and 
of £(x,-%)*, E (x,-) (y,-¥) , and E(y,-y)? for the'fitted 
means’ case are derived. After deriving, in section 2, 
approximate expressions for the determinants involved in 
the generating functions, a general method, due to McGregor 
(1962), for obtaining the joint densities of the regression 
Peeer olenna is presented in section 3. The method is 
carried out in section 4. In sections 5 and 6, approximate 
expressions for the marginal densities and their variances 


are derived. 
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where 2, is the covariance matrix of x. 
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The joint moment generating function of ce Dt, and E* is 
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Die. Evaluation of the Determinants 


_Because the determinants, |A| and la*|, are difficult 
to evaluate exactly, an asymptotic result, first derived by 
McGregor (1959), is used. As the general result is well 
described by McGregor (1959) and by Bielenstein (1963), the 
method is demonstrated here for the special case of approxi- 
mating |A|. Later, |A’| will be approximated by reducing it 
to a scalar multiplied by |A|, 
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Furthermore, Q,7Q5" and Q, are symmetric,and products of 


them are also symmetric and hence commutative. 


Therefore, by (2.2.1) and (2.2.2), 
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Making the substitutions, a, = ~($,+45), and a. = o¢5,(2.2.26) 
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maar} <i Bee 
[(b-2p,) (d-2p,) -U7]fa - A, (a-p,) + A,U] 


[ (b-2p,) (d-2p) - U*] - a(d-2p,) (a-p,) + av 
a(d-2p.,) [b-2p,-atp,] 


ap, (p,-1) (d-29,) ’ 


= [(b-2p,) (d-2p5) - U7} {A,U - A, (e-p5)] 
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YjU[b-2p,-ato,] 


[(b-20,) (4-205) - U1 4U - 234 (c-P5) + yo) 
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To obtain more zeros in the first and last columns of C,, 


in equation (2.2.39), a similar technique is used. 


Let Wy and 5 satisfy the equations, 


and 


Then 


Furthermore, 
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Then 
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Rinall yv,, Ce is multiplied on the right by 
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where [(b-29,) (d-20,) ~ u*)7+ has been factored from the first 


and last columns, and where 


w, = [(b-2p,) (d-295) - U"}lys - uy (2¥g) - uy (2y5)] 
= yg{(b-2p,) (d-205) - U*] + 21 (d-2p,)¥_g + Uva] , 
wy = [(b-2p,) (€-205) - U7} [-1 - uy (a-p,) + QU] 
= -[(b-2p,) (a-295) - U7] + (d-2p,) (a-p,) - i 
= -(d-2p,) [b-2p,-atp,] 
= —9,(p,-1) (d-295)  , 
Ww = [(b-2p,) (d-2p5) - U7} [uyU - uy (e-p,)] 
= -(d-2p,)U + (c-p5)U 
= -Polpo-1)U_ , 
Cha [(b-2p,) (d-2p,) - UI lg - uy (2y, 9) = Uy (2y,,)] 
= yol(b-2p,) (4-299) - U*] + 21(a-20,)¥19 + UY) + 
we = [(b-2p,) (d-2p5) - U"} [-u4(2vg) - u3(2yq) + yg] 
= 2(Uy, + (b-20,)¥5] + vgl(b-2,) (d-2p,) - U] 
We = ((b-2p,) (d-2p5) - U*][-u4(a-p,) + 440] 
= U(a-p,) - (b-29,)U 


= —P, (9, -1)U ; 





= AGP 
wo = [(b-29,) (d-29,) - U"} [u,U ~ Hale-p5) - 1) 
= aif <8 (iO. iean) So Tera (Gere el 
it 2 a os) 
= ~(b-29,) [d-2p-ctp.5] 
= —P5 (95-1) (b-2p,) , 
and 


Wa = [(b-2p)) (d-2p,) - U7 {-uy(2yq9) - U3 (2¥,4) + ¥q9] 


2[Uyy9 + (b-2p,)¥qq] + ¥z9f (b-20,) (d-2p5) - U7] 


The determinant in (2.2.40) may be evaluated asymptotically 
by observing the following: 


Using equations (2.1.19), 

Ye = 49, (p,-1) (d-295) = O(1/n) , 

Y7 = O05(p5-1)U = OU) aay 

Yio™ Y¥201(9,71)U = (B-y°/a) p, (py-1)U SOY Se) 


(8-y7/a) 05 (5-2) (b-2p,) = O(1/n) , 


Yu17 YoPo(Po-1) (b-20)) 


U7} + 21(d-2p,)y, - Uyg] 


W, = Y_l(b-2p,) (d-2p,) 
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+ na(d-20,) [(b-2p,) (d-2p) - U"] 
+ O(1/n) 


= (47S) (4TS-U7)~*[(b-29,) (a-295) - u2}2 


+ 2(d-2p,) [(b-2p,) (d-2p5) - U7] + O(1/n) , (2.2.41) 
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+0(1/n) 


(4rs-u2) (27) 7? 


UL (b-2p,) (a-2p,) - U7) 
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2,2 
O@1/n)) + y,[(b-2p,) (d-2p,) — Ue) 
+ naU[(b-2p,) (a-295) - U7] 
(-27U) (47S-U*) ~*[ (b-294) (d-2py) - U"]? 
+ 2TUL(b-2p,) (d-295) - Cope OUlyni yn) (2.2) 4B) 


-P, (p,-1)0 ’ 


2[Vyy9 + (b-204)¥qq) + yyy l(b-2p,) (4-209) - 071 
O(1/n) + ny, (b-20,) [(b-2p,) (d-295) - U4] 

+ [(b-2p,) (d-295) - 077 

heat 2 
(4rs-v?) (27) ~)(b-2p,) [(b-2p,) (a-2p) - U7] 


+ [(b-2p,) (d-2p4) - uel? + O(1/n) . (2.2.44) 
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If we denote the determinant in equation (2.2.40) by By, 


then JA | = [(b-29,) (d-20,) - u7}"* B, . (2.2.45) 


Ignoring terms which are O(l/n), By may be expanded to the 


2n+1l 
Ww 


Simple form, By av w)B. + (-1) 533 , 
that is By © WB, - weB, A (Qe 240) 
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Bo and B. may again be expanded. Ignoring terms of order oo. 


B. we [Al ; | (2.2547) 


and B, w,lAl : (227A 3) 


where A is the matrix encountered in the *known means' case. 
Hence, using equations (2.2.45) through (2.2.48), 
ja” | = [(b-2p,) (d-2p.) - u7]74 x [wiwe-wew,] |A| (2.2.49) 
= il 2 LenS tet aad ~~! ° yer 
Using relations (2.2.41) through (2.2.44), we may write 
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[w)Wo-wew,] = { 47s] [4ts-v2] 71 [ (b-20,) (d-2,) -U 


+ 20(d-2p,) [(b-2p}) (a-205) - UI} 


(continued) 


ox 
"Rae B to to em + pabtonpl  lebueGes ed whens vee Pie bas of 
e ; t. 7 : 
ce ee 
\ Vee S 
(gt ¢) : |i gw vi pet 






.2B6o ‘aneem mais * Sit af bewssnvaons Kieg0n ate ei a omer. 7 







a 





(8b.5,8) spwordd (EBLE. S) sapanaves “- seca 
Gieyeus) . jal ean «Nie es aaeans tak 









ae 
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+ [4TS-U“] (a-2p,) (b-29,) + 27 (d-2p,) 
2 2 2 
x[ (b-2p,) (d-2p,) = 17) ap 48 [(b-2p,) (d-2p,) =a me 
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and, using equations (2.1.19), 
W.W,-wW.-W, = [(b-29,) (d-29.,) - u2}3 
LB ee a. a al Z 
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where the error in the approximaticn is O(1/n) . 
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Pa BI A Method of Inversion 

The following is a description of a technique for deter- 
mining the joint density function of two ratios with common 
numerator. 

Suppose the ratios are by = E/C and DS = E/D, where C and 
D are positive. The method used involves Manipulating the 
joint moment generating function, 

M(T,S,U) = E[exp(TC+SD+UE) ]. 

The transformation, by =e /C, bo = E/D, u =C+D is made, 
yielding C = bju/ (by tbo) + D = b,u/ (bi +b), be= b bou/(b,+bo), 
with Jacobian, b,b,u2/(b,+b9) >. 
Thus, if the joint density of C, D, Shere) G9 SUEY Sele bi giay bs 


then that of by, bo, and u is 
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and the joint density of by and bo is 
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But M(T,S,U) is the Fourier transform of £(C,D,E) #eand2by 


the Fourier inversion theorem, 


£(C) Uy ee= (2)? {f{Jm(z,s,U) exp[-(TC+SD+UE) ]dudSdT . 
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2.4 The Joint Density of the Two Regression Coefficients 
(a) Known Means 
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Clearly, ~(0) = 0, and the only non-zero term in wy"(0) is that 


corresponding to differentiating the factor bi twice. 
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To find the variance of bi, note that f (b,) was 


obtained from £(b,) by multiplying by expression (2.5.24) 
and changing n to (n-1) . The latter change has no effect 
on the asymptotic variance; the former multiplies (0) by 
expression (2.5.24), and hence, multiplies Var (b,) by the 
Same quantity. 
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CHARTERS ELL 
COMPUTATIONS AND CONCLUSIONS 


3.1 The Calculation and Renormalization of Some Approximate 


Densities with Known Means 


Using equations (2.5.20) through (2.5.23), noting the 
symmetry, one hundred ordinates of each of several densities 
were calculated. Since £(b,) was found to be negligible for 
all values of |b, | larger than 3, these ordinates were calcu- 
lated for b, varying from -3 to +3. The integral, ft (b,)db,, 
was then calculated numerically with a Simpson's rule routine. 
This quantity was printed out, and then used to renormalize 
f(b). Renormalized graphs of £(b,) were then drawn using the 
"Autoplotter" package of the Computing Center. 

A few values of the renormalizing constant have been dis- 
played in Table 1. They exemplify the general trend for the 
constant to approach 1 with increasing n. The value of the 
constant for 07a i> = 0 and n = 30 indicates the accuracy of 
the numerical integration. 

Figures l(a) to (d) illustrate the behaviour of the density 
function for a sample size of 30, and for a representative 
selection of values of the auto-correlations, Py and Po- This 
behaviour conforms Bo what would be expected by considering 


the types of scatter diagrams that are likely to occur. It 
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should also be noted that the density is unchanged if the 


signs of Py and OP, are changed simultaneously, that the vari- 
ance is largest for P1P> large and positive, and that the 
variance is smallest for P1P5 large and negative. All three 

of these tendencies agree with the expression, (2.6.15), for 
the variance of b- 


Figure l(e) illustrates the predicted behaviour of the 


density for increasing values of n. 





Table 1 : Renormalization Constants (RC) : Known Means Case 
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3.2 The Calculation and Renormalization of Some Approximate 


Densities with Fitted Means 


Using equations) (2.5.25) to (2.5.28), once again noting 
the symmetry, one hundred ordinates were calculated and renor- 


* 
malized £or b. varying from —-3 to +3. 


al 

Table 2 contains renormalizing constants corresponding to 
those in Table 1. The same trend for the constant to approach 
l with increasing n is apparent. 

The same pattern in the behaviour of £" (b;) was observed 
as in that of E(b,). Only one set of densities has been 
plotted. PP has been held fixed, and Py allowed to increase 
from 1/3 to 2/3 in steps of 1/12. As would be expected from 

* 


equation (2.6.16) for the variance of bi, the variance de- 


creases monotonically. 
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Table 2 : Renormalization Constants : Fitted Means Case 
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iS The Results of a Monte Carlo Simulation 


To test the accuracy of the approximation, it was decided 
that, rather than run a large and expensive simulation to 
estimate the entire density function, the accuracy of the 
0.025 percentage points should be tested. It is these values 
that would be used in the standard tests and confidence in- 
tervals. 


* 


With Cake =10.5,) four, thousand values of by were gener- 


"D 
acLecdEtoreeach OoLtmenescases, n= 10> 20,eand 30. The propor— 

a © 
CrongoG b's, in each case, greater than the approximate .025 


percentage point was computed. The results, together with the 


statistic used for testing p = 0.025, are presented below. 
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0.02475 


0.02250 





Table 3 : Results of the Simulation 
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Thus, it appears that, for the purpose of calculating per- 
centage points commonly used in testing hypotheses and in con- 
Structing confidence intervals, the approximation is reasonably 
good for samples of size 30 (and perhaps smaller). Larger simu- 
lations considering different values of Ny, Py, and Po and dif- 


ferent percentage points are needed to make this statement 


more precise. 
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Shey (st Conclusions 


The approximate density of the regression coefficient has 
been found to be symmetric about 0, and approximately bell- 
Shaped. The variance was found to depend on Py and Por but 
still to be of order (l/n). Hence, the unbiasedness and con- 
sistency of the least squares estimator is confirmed. Equa- 
tions (2.6.15) and (2.6.16), however, demonstrate that the 
level of the classical tests on Bie and the width of the 


classical confidence intervals are likely to be inaccurate. 
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APPENDIX 
DETAILS OF THE COMPUTATIONS 


Programs were written in Fortran, for use in the Fortran G 
compiler. Double precision variables were found to be needed 
in the calculation of the approximate densities for lo, | or 
|o,| less than 0.1. The calculations are now accurate for los | 
and |o.| as small as 0.01. Graphs were plotted using "Auto- 
plotter. 

The "gamma" function, needed in the calculation of the 
approximate densities when OMe dy 0, was computed using the 
subroutine, "GMMMA", in the "Scientific Subroutine Package". 
Normal random deviates were generated using the subroutine, 
"CS:003A", of the Computing Centre of the University of Al- 
berta. 

Annotated programs, in the order in which their results 


appear in the text, are reproduced on the following pages. 
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C UP TO FIVE APPROXIMATE DENSITY FUNCTIONS ARE 
C CALCULATED (FOR KNOWN MEANS CASE) AND 
C RENORMALIZED. INTEGRATION IS DONE USING 
C SIMPSON'S RULE. THE RENORMALIZING CONSTANTS 
C ARE PRINTED OUT, FOLLOWED BY A MATRIX OF 
C 6 COLUMNS, THE FIRST OF WHICH CONTAINS 
C ABSCISSAE.THE REMAINING COLUMNS CONTAIN 
C CORRESPONDING ORDINATES OF RENORMALIZED 
C DENSITY FUNCTIONS. 
DIMENSION X(100),Y(100,5),N(5),RO1(5),RO2(5),RENOR(5) 
ExMERNALYV,Z, XI 
DIMENSION YT(100) 
C GMMMA 1S THE SSP SUBPROGRAM WHICH EVALUATES 
C THE GAMMA FUNCTION. 
EXTERNAL GMMMA 
C MMAX IS THE NUMBER OF DENSITY FUNCTIONS EVALUATED. 
MMAX=1 
DO 9 1=1,MMAX 
C NesliSaiit “SAMPLE SIZE: 
C ROL 1S THE AUTOCORRELATION OF THE X'S, 
C RO2Z 1S THE AUTOCORRELATION OF THE Y'S., 


DBR EADS ELOOT Nal eo ROI), ROZC)) 
LOS ARORMA Tal Sa2:F 5a 025) 
C K EQUALS HALF THE NUMBER OF ABSCISSAE. 
K=50 
C THE RANGE OF THE ABSCISSAE IS (-OMAX,+0OMAX). 
OMAX=3 
DOs laa = 07 K 
J=50+1 
L=50-| 
Al=l 
AK=K 
C THE X'S ARE DEFINED TO BE EVENLY SPACED OVER (-OMAX,+0MAX). 
X(J)=AlL*OMAX/AK 
X(L)=-X(J) 
DO 1 M=1,MMAX 
G ORDINATES ARE CALCULATED USING THE EVEN SYMMETRY, 
¥(U,M)=F(X(J),RO1(M),RO2(M),V,Z,X1,N(M) ) 
Io aV.Gh, Mavis) 
DO 7 M=1,MMAX 
DOF oT baie. 00 
Tos Weslo aveC lei ) 
C RENOR IS THE AREA UNDER THE DENSITY FUNCTION,F. 
7 RENOR(M)=SIMPCOMAX,X,YT,K) 
C THE RENORMALIZATION CONSTANTS ARE PRINTED OUT. 
Ville C5, 102s) RENOR 
DOF 25 al =05 K . 
DO 2 M=1,MMAX 
J=50+1 
L=50-] 
C THE ORDINATES ARE RENORMALIZED BY 
C DIVIDING BY RENOR, 
¥(JU,M)=Y¥CdU,M)/RENOR(M) 
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2 NGS MEY CS,M) 
HOrS Wali 00 
THE MATRIX DESCRIBED ABOVE IS PRINTED OUT. 
5 WRITES (6,101) BG Cle Glee CL es BYCIL 4) YC, 5) 
101 FORMAT (6F13.6) 
2 hOP 
END 
FUNCTION XI (OMEGA, RO1,RO2) 
SERSEQUAT.I ONS 245 -23- 
Xl =OMEGA**2+SORT( (COMEGA**2+((1-RO1)/(1-RO2) ) **2) 
1* (OMEGA**2+((1+RO1)/(1+RO2) )%*2)) 
RETURN 
END 
FUNCTION VCOMEGA,RO1,RO2,X1) 
SHEREOUATLUON DS? 5.522. 
X1V=X1I (OMEGA, ROI, RO2) 
VE(XIV*(1.-RO2) **24+(1.-RO1) **2) "(XIV 
Le GHEEROZ) ee 2 @VAt RO )**2)/8./RO1/R027 ( 
2X1 V-OMEGA**2) 
RETURN 
END 
BUNGHWON TZ COMEGA,T ROI RO2 XIV) 
SEREREOUATTON= 225221. 
Y=VCOMEGA, RO1,RO2,X1) 
DECROURD 2)" 373; 5 
5 ZElLit+Y+(Y**24+2,*Y) e%,5 
RETURN 
5B L221. +YRr (Yee 242, *Y) x%,5 
RETURN 
END 
BUNCH ONS (OMEGATROIDRO2, V,Z2.X15N) 
SPEMEOUATIONG 225.70. 
PBCROBPSEOSOSUmeeeRO2eeU 00 )eGO a OeS 
XI F=X1I (OMEGA, RO1,RO2) 
VF=V(OMEGA,RO1,RO2,X1I) 
LE= LCUOMEGAS ROL ROZ. xX lev) 
FN=N 
F=38.*FN*(FN-2.)*SQRT(C CABS (RO1*RO2) ) *(1-RO1**2)*(1-RO2**2) 
1)/(3.1415926*(FN-3.))**0.5/(FEN-4&. )/(1.-RO2**2) 
F=Fu( (VF **2+2,*VEF)*(OMEGA**2+((1-RO1)/(1-RO2))**2) 
1*(OMEGA**2+((1+R01)/(1.+R02) )**2))**0.25 
F=F*(1-ZFex2)*(1-ZF)**2/(1-RO1*RO2*ZEF) 
F=F/(X1F*(1.-RO2**2)+(1-ROL**2))/(XIF*(14+RO2**2) 
1+(1+RO1L**2))/CXIF**2*(1-RO2**2) e*2+(1-ROL*4*2) *%2 
2+2.*XIF*((ROL-RO2) **2+(1.-ROL*RO2Z) **2))**0.5 
F=F*(ZF/RO1L/RO2) **((FN-3)/2) 
RETURN 
5 LDECROMCEO. 0.008 GO.0) 4 
RETURN 
4 eCMae. Sale aer api reer Thee as 
RETURN 
6 TN=N+#l1 
TN=TN/2. 
BN=N 
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BN=BN/2. 

CALL GMMMACTN,GT, 0) 

CALL GMMMA(CBN,GB, 0) 
F=(GT/GB/3.1415926**.5)/(1+OMEGA**2)**TN 
RETURN 

END 

FUNCTION SIMPCOMAX,X,Y,K) 

SIMPSON'S RULE NUMERICAL INTEGRATION ROUTINE. 
DIMEN STONSACLOO DY C100) 

SK=K 

H=OMAX/SK 

SUMEND=0. 

SUMMID=0. 

KI=2*K-1 

DOP oma 0, Kine 

L=J 

M=J+1 

SUMEND=SUMEND+Y(L) 

SUMMID=SUMMID+Y (M) 

M=2*K 

SEMP=€2. *SUMEND+4, *SUMMID-YC(0)+Y(M)) *H/3. 
RETURN 

END 
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OD Gia @.@ OicoG. 


UP TO FIVE APPROXIMATE DENSITY FUNCTIONS ARE 
CALCULATED (FOR FITTED MEANS CASE) AND 
RENORMALIZED. INTEGRATION IS DONE USING 
SIMPSON'S RULE. THE RENORMALIZING CONSTANTS 
ARE PRINTED OUT, FOLLOWED BY A MATRIX OF 
6 COLUMNS, THE FIRST OF WHICH CONTAINS 
ABSCISSAE.THE REMAINING COLUMNS CONTAIN 
CORRESPONDING ORDINATES OF RENORMALIZED 
DENSITY FUNCTIONS. 

DIMENSION X(100),¥(100,5),N(5),R01(5), RO2(5),RENOR(5) 
EXTERNAL V,Z,X! 
DIMENSION YT(100) 
CMMMA 1S THE SSP SUBPROCPAM WHICH EVALUATES 
& THE GAMMA FUNCTION. 
EXTERNAL GMMMA 

e MMAX IS THE NUMBER OF DENSITY FUNCTIONS EVALUATED. 
MMAX=5 
pO 9 t=1,MMAX 

g MOUS THE SAMPLE SIZE. 

© RO1 1S THE AUTOCORRELATION OF THE X'S. 

C RO2 1S THE AUTOCORRELATION OF THE Y'S. 


9 READEILOS? NGI )S ROTI), RO2 (1) 
109 PORMATERA SHE 54929 
¢ K EQUALS HALF THE NUMBER OF ABSCISSAE. 
K=59 
C THE RANGE OF THE ABSCISSAE IS (-OMAX,+OMAX). 
OMAX=3 


AK=K 
¢ THce XP GUAREA DEFINED? TO BE EVENLY SPACED OVER (-OMAX,+OMAX). 
X( J) =AI*OMAX/AK 
X(L)=-XCd) 
DO 1 M=1,MMAX 
C ORD. WNATESG ARE, CAUICULIATIDY USI NGHTHE) EVEN SYMMETRY. 
Y(J,M)=F(X(d),ROL(M), RO2(M),V,Z,X1,NCM)) 
¥(L,M)=Y¥ (JU,M) 
DO 7 M=1,MMAX 
NO 13 L=1,1900 
13. FORGE, M) 
C RENOR 1S THE AREA UNDER THE DENSITY EUNCTHOND IPs. 
7  RENOR(M)=SIMPCOMAX,X,YT,K) 
C THE RENORMALI ZATION CONSTANTS ARE PRINTED. OUI 
WRITE (5,101) RENOR 
DO" 92) Mi=O0PiK 
DO 2 M=1,MMAX 
J=50+1 
L=50-1 
C THE ORDINATES ARE RENORMALIZED BY 
C DIN EDANG? BY, (RENOR., 
Y¥(J,M)=Y(J,M)/RENOR(M) 
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YCL,M)=Y(J,M) 

DO 3 1=38,62 

THE MATRIX DESCRIBED ABOVE IS PRINTED OUT. 

NeW Go LOAN) YCIh 1) eY Cl, 2).¥(1,3),YC1,8),Y C155) 
FORMAT (6F13.6) 

STOP 

END 

FUNCTION XI (OMEGA, RO1,RO2) 

SEE EQUATION 2.5.23. 

X 1 =OMECA**2+SORT ( (OMEGA**2+((1-RO1)/(1-RO2)) **2) 
1* (OMEGA**2+((1+RO01)/(1+RO2))**2)) 

RETURN 

END 

FUNCTION VCOMFGA,RO1,R02,X1) 

SFE EQUATION 2.5.22. 

XIV=X1COMFGA, ROL, RO2) 
V=(XEV*(1.-RO2)*#24(1.-ROL)**#2)* (XIV 
1*(1.+RO2)**24+(1.+RO1)**2)/8./RO1/RO2/ ( 
2X1V-OMEGA**2) 

RETURN 

END 

FUNCTION ZCOMEGA, RO1,RO2,X1,V) 

SEE EQUATION 2.5.21. 

Y=V (OMEGA, RO1,RO2,X1) 

IF(RO1*RO2) 3,3,5 

Z=1L.tY+(V¥ee242 4Y) eS 

RETURN 

Z=1L.+Y¥~(Y¥e*24+2,%Y) #5 

RETURN 

END 

FUNCTION FCOMEGA,RO1,RO2,V,Z,XI,N) 

SEE EQUATION 2.5.20. 
hacRoimeav020n0R. PO? 60.060) GO 1055 
XIF=X1(OMFCA, ROI, RO2) 

VF=V (OMEGA, RO1,R02,X1) 

ZF=Z(OMEGA, RO1,RO2,X1,V) 

FN=N-1 
F=8.*FN*(FEN-2,)*SORT( (ABS (ROI*RO2)) *(1-RO1L¥*2)*(1-RO02**2) 
1)/(3.1415926*(FN=3.))**0.5/(FN-b.)/(1.-RO2**2) 
Feb ((VFee242,*VE)* (OMECA®*2+((1-P01)/(1-RO2)) **2) 
L* (OMEGA**2+((1+RO1)/(1.+RO2)) **2))**0,25 
Fafe(1-ZF**2)*(1-ZF)**2/(1-RO1*RO2*ZF) 
F=E/(XIE*(1,-RO2**2)4+(1-ROL#*2))/(XIF*(1+RO2**2) 
1+(1+ROL**2))/ (XI Fe*2*(1-ROQ*H2 )H*24(1-ROL ERD) 42 
24+2.*X1F*((ROI-RO2Z)**24+(1.-ROL*RO2) **2))**0.5 
F=F*(ZF/ROL/RO2)**((FN-3)/2) 
F=F*(1+ZF)/(1-RO1)/(1-RO2) 
F=E*(XIF*(1-RO2)**2+(1-RO1) **2) 
F=E/(XIF*(1+RO2**2)+(1+RO1**2)) 

RETURN 

PoueOtmED nO. 0M G0.41.0. 4 

RETURN 

1E(RO2.E0.0.0) GO TO 6 

RETURN 
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TN=N 

TN=TN/2. 

BN=N-1 

BN=8N/2,. 

CALL GMMMACTN, GT, 0) 

CALL GMMMACBN, GB,0) 

F=(GT/GB/3.1415926** .5)/C1+OMEGA**2 )**TN 
RETURN 

END 

FUNCTION SIMPCOMAX,X,Y,K) 

SIMPSON'S RULE NUMERICAL INTEGRATION ROUTINE. 
DEMENS TON X GED0De Y (100) 

SK=K 
H=OMAX/SK 
SUMEND=0. 
SUMMID=0. 
K1=2*K-1 

DOm8s J=0,K1,2 
L=J 

M=J+1 
SUMEND=SUMEND+Y(L) 

SUMMID=SUMMID+Y(M) 

M=2*K 
SLMP=(2.*SUMEND+4,*SUMMID-Y(0)+Y(M))*H/3. 
RETURN 

END 
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DIMENSION 8(4000) 
EATERNAL €5003D 
Peon ser) TXT) NS 


Cc XT 1S AN ARBITRARY INTEGER REQUIRED BY CS003A 
C TO START GENERATING NORMAL DEN Ie Si. 
C NS JS THE NUMBER OF REGRESS 101 COEETIGIENTS CALCULATED. 
CALE CSOOSACTXT) 
BreCS 7107) NZ ROIFRO2,8C 
C N IS THE Serie MeLeasey TO CALCULATE FA REGRESSION COEFFICIENT, 
C ROL AND RO2 ARE THE py LO GORBEL ATLONS:. 
G POmlomine .025 PERCENTAGE POINT OF THE APPROXIMATE DENSITY FUNCTION 
D7 eORMAT C12, OFS, 27 Lo. 6) 
101 FORMAT lis LS) 
DO 1 K=1,NS 
1 BCK)=BMC(N,RO1,R02) 
& BelS Ay VECTOR OF LENGTH NS CONSISTING OF 
C THE GENERATED REGRESSION COEPEI GIENTS. 
BR=0, 
O05 1=),NS 
[ECBC!).GT.BC) BR=eRR+1, 
C BR 1S THE NUMBER OF REGRESSION COEFFICIENTS GREATER THAN 8c. 
8 CONTINUE 
FNS=NS 
BPR=BR/FNS 
C mom Us THE PROPORTION OF REGRESSION COEFFICIENTS GREATER THAN Bc. 


PATI ECG,108) BPR 
108 FORMAT(F1N.6) 

3s hOP 

END 

FUNCTION BMC(N,RO1,R02) 

THE FUNCTION BMC GENERATES REGRESSION COEFELGEENTS 

PebeohMPEES OF SIZE UP TO 50. 

DIMENSION X(50),Y(50) 

DOe2 het) hl 

WalieCs0050(0.,1.,U,V,KN) 

amet) COTO 3 

J=1-1 

XCE)=RO1L*X(J) +U 

YC1)=RO2*Y( J) +V 

mre OL) YL) 

XX=XX+X(1) 

etme yey. C18) 

XXX=XXX+X(C])**2 

GOl 1022 
3 X(1) =U 

Y(1)=V 

XY=XCI 


Go) 


XXX=X (CJ )**2 
2 CONT] NUE 

BMC=(N*XY- “HX*YY) / ON*XX yoy yun) 
102 FORMAT (F13.6) 
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